%Define a function which returns the residual between your matrix and your fitted curve
ms=123:127;
fs=0.20:0.03:0.8;

x=kron(ms,ones(21,1));
y=kron(fs,ones(5,1))';
myfun = @(params) params(1)+ 1/(2*(1-params(2)^2))*(((x-params(3))/params(4)).^2+((y-params(5))/params(6)).^2-2*params(2).*((x-params(3))/params(4)).*((y-params(5))/params(6))) - resultsD000R1;

%Define initial guesses for parameters a, b
params0 = [360.6,1296,0.5425,1198,0.1817,0.9307];
%Add lots of debugging info
opts = optimset('Display','Iter');
%Fit
fitparams = lsqnonlin(myfun,params0,[],[],opts);

params0=fitparams;
myfun2 = @(z) params0(1)+ 1/(2*(1-params0(2)^2))*(((z(1)-params0(3))/params0(4)).^2+((z(2)-params0(5))/params0(6)).^2-2*params0(2).*((z(1)-params0(3))/params0(4)).*((z(2)-params0(5))/params0(6)));

%optimoptions(@fminunc,'TolFun',1e-100);
[min, minval] = fminunc(myfun2,[160,0.1])